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We investigate the Green function of two-dimensional dense random packings of grains in order 
to discriminate between the different theories of stress transmission in granular materials. Our 
computer simulations allow for a detailed quantitative investigation of the dynamics which is difficult 
to obtain experimentally. We show that both hyperbolic and parabolic models of stress transmission 
fail to predict the correct stress distribution in the studied region of the parameters space. We 
demonstrate that the compressional and shear components of the stress compare very well with the 
predictions of isotropic elasticity for a wide range of pressures and porosities and for both frictional 
and frictionless packings. However, the states used in this study do not include the critical isostatic 
point for frictional particles, so that our results do not preclude the fact that corrections to elasticity 
may appear at the critical point of jamming, or for other sample preparation protocols, as discussed 
in the main text. We show that the agreement holds in the bulk of the packings as well as at the 
boundaries and we validate the linear dependence of the stress profile width with depth. 

PACS numbers: 81.05.Rm, 81.40.Jj 



I. INTRODUCTION 

The theoretical understanding of the structural and 
mechanical properties of static granular assemblies is still 
an open issue 1]. Indeed, there is no consensus on how 
the stress distribution should be expressed and how gran- 
ular materials respond to applied perturbations. Exper- 
imental measurements of stress distributions in silos and 
sandpiles 0, 0, 0, W\ seem to be incompatible with a 
simple elastic theory and have led to the development of 
new mechanical formulations for granular matter where 
new assumptions are proposed |El3- 

The equations of equilibrium for the stress tensor, cr^j. 
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where p is the density of the medium and g the gravi- 
tational acceleration vector. From this equation we see 
that the stress is indeterminate. For instance in 2D, it 
provides two independent equations for the three compo- 
nents of the stress tensor (the stress tensor is symmet- 
rical, (Jxy = cFyx)' Thus a third constitutive equation 
(the missing equation) characterizing the behavior of the 
material is needed in order to provide closure to the prob- 
lem. In the elastic formulation closure is provided by the 
introduction of the strain field and new constitutive re- 
lations relating the strain with the strain (Hooke's law). 
This leads to elliptic equations governing the stress dis- 
tribution Q. 

Unlike elasticity — where the closure is provided by 
Hooke's law — new granular models propose closure rela- 
tions involving only stress tensor components . Hyper- 
bolic equations of the same mathematical structure as a 
wave equation are thus obtained. As a consequence, a 
load applied at the surface propagates along two charac- 
teristic directions. 



Another theoretical framework of stress transmission 
in granular matter (the q-model) was introduced by Cop- 
persmith et at. 7]. It considers the transmission of in- 
terparticle forces by a random redistribution of forces 
among nearest neighbors. At the continuum scale j^, 
the stress components satisfy a parabolic equation of dif- 
fusion, expected from the stochastic uncorrelated nature 
of the transmission process. 

Thus, different theories predict different stress distri- 
butions according to the nature of the equations govern- 
ing the stress. The importance of resolving which of these 
models apply under which circumstances is of major the- 
oretical and practical interest in order to reveal the real 
mechanical properties of static granular assemblies. 

In order to test the validity of the different models, a 
simple experiment was proposed by de Gennes : the 
calculation of the Green function by applying a pertur- 
bation force on a grain in the granular medium and mea- 
suring the resulting stress field perturbation. The elastic 
and the new models predictions are then confronted. In 
the most simple case of an homogeneous isotropic ma- 
terial, the elastic framework predicts that the vertical 
stress profile resulting from the perturbation is a single 
peak centered below the perturbation with a width W 
proportional to the depth y from the perturbation. In 
the hyperbolic models, the stress profile presents a dis- 
tinctive double peak and in the parabolic models the pro- 
file is a single peak as in the elastic prediction but with 
a width proportional to ~ y^. 

The Green function problem has divided the granular 
community into two. For some researchers, the prob- 
lem is well settled by now, in particular after the exper- 
imental work of Reydellet and Clement [n| and Ceng 
et al. p^] in favor of the validity of elasticity theory 
to model the response function. Their experiments in- 
dicate a single peak in the stress distribution in cjyy in 
qualitative agreement with the elastic framework. More- 
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over, 1 shows that the elastic framework can be further 
extended to anisotropic system. However, other experi- 
ments 111 111 111 El El, El 113 show evidence for the 
possible validity of alternative approaches: 14] discusses 
the lateral diffusion of vertical stress increments due to 
the local application of a force in a brick wall; El HI 
focus on relatively small scale response experiments and 
^201 on ordered lattices; The study on isostatic networks 
of 13 and 13 show deviations from elasticity. These 
groups argue that on approaching the isostatic point the 
elasticity theory would break down and new approaches 
need to be considered to describe the stress response. 
On the other hand some groups |21j argue that in the 
isostatic limit the hyperbolic equations do not describe 
the stress field. Recent work ^ suggest that there ex- 
ist a characteristic length above which elasticity is valid 
and below anomalous behaviour of isostaticity should be 
observed. Also, it was shown that the degree of ordering 
23, 24, 2I and the texture properties of packings 
|28j affect strongly the stress transmission. Even though, 
there is general consensus that the elastic theory may de- 
scribe qualitatively the system, many discrepancies still 
remain |2^. A microscopic study is needed to perform 
a quantitative assessment of the validity of the different 
theoretical approaches. Simulations offer an optimum 
scenario to compare the different theories and allow to 
obtain microscopic information 2£., 29] and the exact 
system dynamics [sl, difficult to observe in experiments. 

Here we perform MD simulations to gain more micro- 
scopic insight into the mechanical response of granular 
matter. Our results indicate that the elastic theory can 
describe the simulations better than the other models. 
Indeed, we find a stress response whose vertical compo- 
nent shows one peak. Furthermore we show that elastic- 
ity predicts very accurately all components of the stress 
tensor at all depths in the packing, for all the ranges of 
pressures, porosities and frictional properties studied in 
this work. 



II. NUMERICAL SIMULATIONS. 

We developed a numerical code based on the Dis- 
crete Element Method (DEM). The model deals with 
an assembly of elastic spheres interacting via the Hertz- 
Mindlin contact laws. The grains interact with one an- 
other via non-linear Hertz normal forces which de- 
pend on the overlap between two spheres, and frictional 
transverse forces which depend both on the shear and 
normal displacements between the grains: 

fc 2c i?V2^3/2 
Af- = Ct{Rwy/^As, 

where R = 2RiR2/{Ri + R2) with Ri and R2 the radii 
of the spheres in contact. The normal overlap between 
the two grain is and As is the relative shear dis- 
placement between the two grain centers. The constants 
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FIG. 1: System geometry. The left part of the figure shows a 
detail of the bidisperse particle packing and the typical grid 
used to calculate the coarse-grained stress profile, while the 
right part shows the heterogeneity in the interparticle forces 
("force chains"). A small force / is applied on a single grain 
located at the center of the sample (origin of reference axes). 
The vertical boundaries are composed of rigid grains while 
the horizontal boundaries are periodic. 



Cn = 4G/ (l — i') and Ct = 4G/ (2 — v) are defined in terms 
of the shear modulus G = 29GPa and the Poisson's ratio 
u = 0.2 of the material of the grains. Coulomb friction, 
ft ^ f^ffn^ .31] is also included in the model and we set 
the friction coefficient /i/ = 0.3. We also include viscous 
damping terms in the equations of motions (translational 
and rotational) to allow the system to relax toward static 
equilibrium. In this study, the grains spherical are con- 
strained to move in 2D. 

The packings are prepared using binary mixtures of 
particles at equal concentrations characterized by a ra- 
dius R = (0.1 ± 0.01) mm in order to avoid crystal- 
lization. We study different sample sizes ranging from 
N = 10000 to N = 250000 particles equilibrated using 
periodic boundary conditions in the horizontal directions 
and two rigid walls made of grains at the top and the bot- 
tom of the system (see Fig. QJ. In order to prepare the 
larger samples {N > 10000) at static equilibrium and 
save computational time on the preparation procedure, 
we repeated periodically in space the elemental sample 
(A^ = 10000) equilibrated using periodic boundary con- 
dition in both horizontal and vertical directions. For in- 
stance, to prepare a packing with size N = 40000, we re- 
peat the elemental packing twice in both horizontal and 
vertical directions. We found no significant difference for 
large system size, so that most of the results presented 
here are for 10000 particles, which were obtained without 
applying the repeating technique. We use a numerical 
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compaction protocol designed to prepare confined dense 
granular isotropic packings explained in detail in 32]. 
We set the gravity g = since the pressures of confine- 
ment considered in this work are sufficiently high. 

A packing is isostatic ^34^ when the number of contact 
forces equals the number of force balance equations. In a 
packing of perfectly smooth (frictionless) particles, there 
are NZ/2 unknown normal forces and DN force balance 
equations, where D is the dimension and Z is the average 
coordination number. This result in a minimal coordina- 
tion number needed for mechanical stability as Zc = 2D, 
i.e., Zc = 4 and 6 in 2D and 3D (See Fig. [2 the system is 
at the isostatic limit with Z :^ Zc as p ^ 0). In the case 
of packings of perfectly rough particles, which is realized 
by frictional particles with infinite friction /i ^ oo [note 
that in this case there are still tangential forces given 
by the Mindlin elastic component Eq. 2], in addition 
to NZ/2 unknown normal forces and DN force balance 
equations, there are {D — l)NZ/2 unknown tangential 
forces and D{D — l)N/2 torque balance equations |35| . 
Thus the coordination number in the isostatic limit is 
Zc = + 1, i.e., Zc = 3 and 4 in 2D and 3D for frictional 
packings. In such isostatic packings, there is possibly 
a unique solution for the forces between particles for a 
given geometrical configuration, because the number of 
equations equals the number of unknowns. The existence 
of isostaticity is the foundation of recent theories of stress 
propagation in granular materials 35, 36, 37]. 

Here we use a specific preparation protocol which sets 
the friction between particles to zero. After a sample 
is equilibrated at a given pressure, we consider two sit- 
uations to calculate the Green function: (i) frictionless 
case, /i = 0; (ii) frictional case, /i = 0.3. Note that 
in (ii) friction is turned on an equilibrated frictionless 
packing only for the calculation of the Green function. 
Thus, this case is far from the isostatic frictional case 
Zc = + 1 38]. We use a frictionless preparation be- 
cause (1) we want to avoid issues of path dependence in 
the sample preparations and (2) to obtain the isotropic 
packings with large coordination number (6 in 3D, see 
|3^, l40|). Our preparation protocol tries to mimic the 
shaking that is usually applied to the packings in exper- 
iments to remove unstable voids that would otherwise 
render the system unjammed. In a sense, the packings 
that we obtain are analogous to the " reversible" packings 
obtained in the compaction experiments of Nowak et at 
j41j. This preparation scheme was tested against exper- 
iments in our previous computations of sound speeds in 
granular materials ^3^ . We find very good agreement be- 
tween experiments and simulations. Thus we believe that 
it is a reliable method to prepare the packing. However 
it will be very interesting to see the response behavior 
for a packing prepared with friction. For this case, ex- 
treme care has to be taken to assure that the packing is 
jammed and reversable, i.e., that there are no unstable 
voids. We have recently developed a preparation pro- 
tocol to provide jammed frictional packings |33- The 
Green function calculation will be applied to this fric- 
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FIG. 2: Coordination number versus pressure obtained in 
2D frictionless packings, which follows the scaling behavior: 
Z{p) - Zc ^ ^o.27±o.o7 ^^^Yi Zc = 3.8. Here the system is 
isostatic with Z ^ 4 as p ^ 0. Note that the critical co- 
ordination number Zc at isostatic limit of frictionless case is 
actually a little smaller than 4, because of the existence of 
floaters which carry no forces. 



tional packings to investigate the validity of elasticity as 
the frictional isostatic point is reached. 

We check the isotropy and randomness of the pack- 
ings by calculating both the texture tensor and the 2D 
orient ational order parameter 42]. In a first part, the 
simulated granular aggregates are prepared with target 
pressure of 2.10^ N/m (equivalent to a pressure of lOMPa 
in 3D) and result with an average coordination num- 
ber Z ^ 4.3 (see Fig. |2) and a solid volume fraction 
(j) ^ 0.860. In a next part, in order to be closer to the 
isostatic point, we show results for lower pressures 2.10^ 
N/m and 2.10^ N/m (equivalent to lOOKPa and IMPa in 
3D) with Z ^ 3.9 and 4.1 (see Fig. |2| compared to the 
isostatic limit Zc = 4, and (j) ^ 0.838 and 0.842 slightly 
above the random close packing limit of (j) ^ 0.82. The 
next step of the preparation protocol consists to break 
the vertical symmetry of the system by removing the pe- 
riodic conditions in the vertical direction and replacing 
it by rigid walls made of rough particles at this time. 

Next, we select the particle nearest to the center of the 
system and apply a vertical downward force. The ampli- 
tude of the force, / = 10~^(/) (where (/) is the average 
normal force in the packing), is small enough to assure 
that we measure the linear response regime. To calculate 
the local stress tensor aij{x, y), the system is subdivided 
into square cells of size / = 0.5 mm containing about 4 
to 5 particles depending on their local arrangements (see 
upper left corner of Fig.Q). The results are independent 
of the size of the coarse-graining cell in the range / to 3L 

We have performed these simulations using three dif- 
ferent system sizes in a range L to 3L in order to find the 
minimal size to consider to avoid finite size effects. The 
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aspect ratio 1 to 3 is needed for the appropriate study of 
the stress response function because the width is of the 
order of y. We find that the vertical response is less sen- 
sitive than the two others component to the sample size. 
We have also compared the vertical response functions 
above and below the perturbation, respectively in the di- 
latation region and in the compression region. We find an 
horizontal symmetry (Jii{x^ y) ^ (Jii{x^ —y)i showing that 
the compression response and the dilatation response are 
almost identical in the studied linear regime. 

We measure (x, y) in units of / with x = 0, = at the 
center of the packing. We consider system sizes ranging 
from L = 40/ (TV = 10000) to L = 120/ (TV = 90000). We 
find no appreciable size effects in this range for (Jyy{x^ y) 
and slight size effects for cFx^i^-, y) and cr^yix^ y) saturat- 
ing above L = 80/ {N = 40000). The perturbation of 
the stress is calculated as aij = afj^ — (7^^^ , where a^J-^ 
is the stress of the initial reference state before the per- 
turbation and crfj^ is the stress after the perturbation is 
applied {cju > corresponds to compression). Because 
of the discrete nature of the material and the strong in- 
homogeneity of the contact network 42] (see Fig.^, we 
average the stress response functions over many configu- 
rations of particles (about 20). 



III. WAVE PROPAGATION AND STRESS 
DISTRIBUTION. 



The application of the force on the central particle gen- 
erates a pressure wave which propagates in the granu- 
lar packing (as seen in the inset of Fig. (Ht) and is re- 
fiected at the rigid and periodic boundaries. The wave 
is dissipated slowly over time until the packing reaches 
a new static state at mechanical equilibrium. Figure OK 
presents the evolution with time of the vertical stress 
(Jyy{x = 1,?/, t > 0) at different depths y below the per- 
turbation in the compression region. At all depths, a 
first stress maximum is observed at the passage of the 
front wave (local instantaneous response in the propaga- 
tive regime, point A), followed by a rapid oscillation and 
then a slow relaxation (point B). Then, a new stress peak 
of smaller amplitude is measured (point A') which corre- 
sponds to the wave refiected at the boundary followed by 
the same decrease sequence towards complete relaxation 
to the new static state (point B'). The further the point 
of measurement is, the smaller is the amplitude of the 
first stress maximum resulting from the expansion of the 
front. 

Figure Ob plots the stresses at equilibrium as well as 
the perturbation of the normal interparticle forces, F^. 
We clearly see the existence of a single peak in a, ,,, in 
agreement with elasticity and the experiments of [ll|, [13 
and in disagreement with the predictions of the hyper- 
bolic theories of stress transmission. 
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FIG. 3: (a) Time evolution of Gyy in a single frictional pack- 
ing at X = and at four different depths y (as indicated by 
the squares in the inset, showing the stress propagation). We 
observe an instantaneous stress peak (point A) in the prop- 
agative regime as a result of the wave front first passage fol- 
lowed by a relaxation to point B to the static regime and a 
subsequent peak followed by another relaxation (from A' to 
B'). The slow damping drives the system from a propagative 
regime to a static state, (b) Typical stress response functions 
of one packing for the three components of the stress tensor 
c^ccx, cFyy and (Txy and the enlarged detail of the perturbation 
of the normal forces Fn on the contact network shown in Fig. 

m 



IV. COMPARISON WITH ELASTICITY. 

Next, we make a detailed comparison with the predic- 
tion of elasticity theory. In the case of an infinite elastic 
medium the solution was calculated by Boussinesq : 



1 



(i + (f)^) 
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(3) 



However, there is no analytical solution for the boundary 
conditions used in the present part of the study, i.e., rigid 
top and bottom rough boundaries where the strain satis- 
fies Uy{x^ y = ±i^/2) = 0, and Ux{x^ y = ±i^/2) = (this 
conditions can be expressed in terms of the stress ten- 
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FIG. 4: (a) Comparison between the solutions for the stress 
between simulations (symbols) and elasticity (lines) at differ- 
ent depths y and at a fixed pressure. The main plot compares 
the vertical component of the stress tensor ayy while the in- 
sets compare the two others components Gxx and Gxy. (b) 
Half-width amplitude of the averaged vertical stress profiles 
with increasing depth. A linear dependence \s W ^ 1-15^ is 
found, corroborating the elastic solution W ^ l.SO^/. 
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FIG. 5: Comparison between the vertical stress prediction 
of elasticity and the simulations at different depths y and 
for five different pressures ranging from p — 2.10^ N/m to 
p — 2.10^ N/m for (a) rough boundaries conditions / fric- 
tional particles, (b) smooth boundaries conditions / friction- 
less particles. The insets show the half- width amplitude of 
the corresponding stress profiles with increasing depth. Lin- 
ear dependence W ^ y is found. 



sor as dyaxx{x, y = ±L/2) = (2 + i^)dxcrxy{x, y = ±L/2) 
and (Jxxix, y = ^L/2) = v(jyy{x^ y = -^L/2) respectively, 
where v is the Poisson ratio of the entire packing [271 ) and 
horizontal periodic boundary conditions. Then, we find 
the solution of Eq. Q numerically 43] . The applied per- 
turbation Q{x) is taken to be a narrow Gaussian centered 
at X = 0, y = 0, such that (Jyy{x^ y = 0, t > 0) = Q{x). 

In the theory, the stress distribution is not sensitive to 
the Poisson ratio as we find numerically. Moreover, this 
result was confirmed in previous simulations Since 
the results are independent of a possible choice is 
to set the Poisson ratio equal to the effective Poisson 
ratio, v = Ve 0.28, which is given by simulations 
[3^ . Moreover, the effective Poisson ratio of a granu- 
lar packing is independent on pressure. [Note that some 
groups found the Poisson coefficients varying with confin- 
ing pressure, for frictionless packings at the approach of 
the isostatic limit 33] and frictional packings with lower 
coordination number 28]]. This can be seen from the 
Effective Medium Theory calculations [s^, which gives: 
Ve = {Ke — 2/3/ie)/(2i^e " 2/3/ie), whcrc Ke and fie are 
the effective bulk modulus and shear modulus, respec- 
tively. Due to the same power law relations of and 
lie versus pressure, ~ p^^^ and /ig ~ p^^^^ the ef- 



fective Poisson ratio Ve is independent of the pressure. 
Therefore we use the same z/g for all our systems at dif- 
ferent pressures. We note, however, that the results of 
[3^ are for 3D while our systems are 2D. Other choices 
of V would give the same results as presented here. Fig- 
ure^ shows how the elastic solution compares with the 
numerical results for the stress components. The elastic 
solution provides a very satisfactory fit at all depth for 
all the components of the stress tensor. 

Width of the Stress Profile. — We obtain the width of 
the vertical stress profile as a function of the depth y by 
calculating the half-amplitude spread of the distribution 
Figure 2b shows a linear dependence on the depth, 
W ^ y/m general agreement with elasticity (and in dis- 
agreement with the parabolic model 7]). The numerical 
response (W = 1.15?/) exhibits a slightly narrower re- 
sponse than the elastic^lution {W = 1.30?/); it was also 
found in experiments [llj and could be attributed to the 
disorder of the packing. 

Pressure and Friction Effects. — We also study the 
stress response functions for different pressures (from 
p = 2.10^ N/m to p = 2.10^ N/m) and porosities on 
smaller systems {L = 80/ and N = 20000) shown in Fig. 
[5Ji. We find that the stress profiles averaged over 20 real- 
izations are in extremely good agreement with the elastic 
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prediction for rough boundaries. Moreover the profiles 
do not show appreciable pressure dependence. This is 
in further agreement with elasticity. Indeed the elastic 
solution depends only on the Poisson ratio v which is 
constant with pressure 32]. Moreover, we find the same 
linear dependence with the depth (W = 1.20?/) for the 
studied range of pressure (see inset of Fig. [5^). The 
agreement with elasticity holds very well also for the two 
others components. 

We also investigate the effect of friction by perform- 
ing simulations with frictionless packings (/ij = 0) in 
the same range of pressure and with the same system 
size shown in Fig. [SJd. This time, we compare the 
stress profiles averaged over 50 realizations (due to much 
larger fiuctuations) with the elastic prediction for smooth 
boundaries (i.e., crxy{x^y = ±i^/2) = 0). We find that 
the stress profiles show a small pressure dependence and 
that the agreement with the elastic prediction holds bet- 
ter for the smallest pressure. The reason for the dis- 
crepancy could come from the fact that even though the 
particles are frictionless, the fixed grains at the top and 
bottom surface induce a roughness, the effect of which is 
expected to be more important at high pressures. Indeed 
despite non-frictional particles, a non-zero shear stress is 
found at the boundary, whose amplitude increases with 
pressure. The linear dependence of the half- amplitude 
show similar dependence with pressure. 



ing into account the exact geometry of the simulations, 
we compare the numerical results and find very good 
agreement with elasticity for the three components of the 
stress tensor in the frictional case. Whereas, experiments 
probe only the vertical stress profile at the boundary, our 
numerical simulations are able to show the agreement 
with elasticity at all depth inside the packings and all 
the components. Furthermore, we confirm that pressure 
is not a critical parameter for the stress response function 
of frictional packings. 

However, it is important to note that we have not in- 
vestigate the frictional isostatic point. Near this critical 
point the system may show deviations from elasticity. 
Thus, we may still expect corrections to appear if the 
system is prepared closer to the isostatic frictional jam- 
ming point Zc = 3 in 2D. It will be of interest, then, to 
continue this work and investigate the critical regime as 
well. We also demonstrate that the response of friction- 
less packings is consistent with elasticity; it exhibits pres- 
sure dependency and tends to agree with elasticity at low 
pressure. We attribute this effect to the roughness of the 
boundaries in our simulations. In further work, it would 
be interesting to confront the more general anisotropic 
elastic theory to the response of loose and anisotropic 
packings obtained by other compaction protocols. 



V. SUMMARY 

Our numerical work shows that, for the studied fric- 
tional and frictionless packings, the vertical component 
of the stress response has a single peak at all depths and 
that the half-amplitude width of the vertical stress pro- 
files is proportional to the depth. By calculating the 
stress profiles using the elasticity framework and tak- 
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